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Summary 

X-ray  computed  microtomography  (XMT)  is  used  for  high- 
resolution,  nondestructive  imaging  and  has  been  applied  success¬ 
fully  to  geologic  media.  Despite  the  potential  of  XMT  to  aid  in 
formation  evaluation,  currently  it  is  used  mostly  as  a  research  tool 
One  factor  preventing  more  widespread  application  of  XMT  tech¬ 
nology  is  limited  accessibility  to  microtomography  bcamlincs.  An¬ 
other  factor  is  that  computational  tools  for  quantitative  image 
analysis  have  not  kept  pace  with  the  imaging  technology  itself 

In  this  paper,  we  present  a  new  grain-based  algorithm  used  for 
network  generation.  The  algorithm  differs  from  other  approaches 
because  it  uses  the  granular  structure  of  the  material  as  a  template 
for  creating  the  pore  network  rather  than  operating  on  the  voxel  set 
directly.  With  this  algorithm,  several  advantages  emerge:  the  al¬ 
gorithm  is  significantly  faster  computationally,  less  dependent  on 
image  resolution,  and  the  network  structure  is  tied  to  the  funda¬ 
mental  granular  structure  of  the  material.  In  this  paper,  we  present 
extensive  validation  of  the  algorithm  using  computer-generated 
packings.  These  analyses  provide  guidance  on  issues  such  as  ac¬ 
curacy  and  voxel  resolution.  The  algorithm  is  applied  to  two  sand¬ 
stone  samples  taken  from  different  facies  of  the  Frontier  Formation 
in  Wyoming,  USA,  and  imaged  using  synchrotron  XMT.  Morpho¬ 
logic  and  flow-modeling  results  are  presented 

Introduction 

Subsurface  transport  processes  such  as  oil  and  gas  production  arc 
multisealc  processes.  The  pore  scale  governs  many  physical  and 
chemical  interactions  and  is  the  appropriate  characteristic  scale  for 
the  fundamental  governing  equations.  The  continuum  scale  is  used 
for  most  core  or  laboratory  scale  measurements  (e.g.,  Darcy  ve¬ 
locity,  phase  saturation,  and  bulk  capillary  pressure).  The  field 
scale  is  the  relevant  scale  for  production  and  reservoir  simulation. 

Multisealc  modeling  strategics  aim  to  address  these  complexi¬ 
ties  by  integrating  the  various  length  seales.  While  pore-scale  mod¬ 
eling  is  an  essential  component  of  multisealc  modeling,  quantita¬ 
tive  methods  arc  not  as  well-developed  as  their  continuum-scale 
counterparts.  Hence,  pore-scale  modeling  represents  a  weak  link  in 
current  multisealc  techniques. 

The  most  fundamental  approach  for  pore-scale  modeling  is 
direct  solution  of  the  equations  of  motion  (along  with  other  rel¬ 
evant  conservation  equations),  which  can  be  performed  using  a 
number  of  numerical  techniques.  The  finite-clement  method  is  the 
most  general  approach  in  terms  of  the  range  of  fluid  and  solid 
mechanics  problems  that  can  be  addressed.  Finite-difference  and 
finite -volume  methods  are  more  widely  used  in  the  computational 
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fluid  dynamics  community  The  boundary  element  method  is  very 
well  suited  for  low-Rcynolds  number  flow  of  Newtonian  fluids  (in¬ 
cluding  multiphase  flows).  Finally,  the  lattice-Boltzmann  method 
has  been  favored  in  the  porous-media  community  because  it  easily 
adapts  to  the  complex  geometries  found  in  natural  materials. 

A  less  rigorous  approach  is  network  modeling,  which  gives  an 
approximate  solution  to  the  governing  equations.  It  requires  dis¬ 
cretization  of  the  pore  space  into  pores  and  pore  throats,  and  trans¬ 
port  is  modeled  by  imposing  conservation  equations  at  the  pore 
scale.  Network  modeling  involves  two  levels  of  approximation. 
The  first  is  the  representation  of  the  complex,  continuous  void 
space  as  discrete  pores  and  throats.  The  second  is  the  approxima¬ 
tion  to  the  fluid  mechanics  when  solving  the  governing  equations 
within  the  networks.  The  positive  tradeoff  for  these  significant 
simplifications  is  the  ability  to  model  transport  over  orders-of- 
magmtude  larger  characteristic  scales  than  is  possible  with  direct 
solutions  of  the  equations  of  motion.  Consequently,  the  two  ap¬ 
proaches  (rigorous  modeling  of  the  conservation  equations  vs.  net¬ 
work  modeling)  have  complementary  roles  in  the  overall  context 
of  multiseale  modeling  Direet  methods  will  remain  essential  for 
studying  first-principles  behavior  and  subpore-scale  processes 
sueh  as  diffusion  boundary  layers  during  surface  reactions,  while 
network  modeling  will  provide  the  best  avenue  for  capturing  larger 
characteristic  scales  (which  is  necessary  for  modeling  the  pore-to- 
eontinuum-sealc  transition). 

This  research  addresses  one  of  the  significant  hurdles  for  quan¬ 
titative  network  modeling:  the  use  of  high-resolution  imaging  of 
real  materials  for  quantitative  flow  modeling.  We  focus  in  particu¬ 
lar  on  XMT  to  obtain  3D  pore-scale  images,  and  present  a  new 
technique  for  direct  mapping  of  the  XMT  data  onto  networks  for 
quantitative  modeling.  This  direet  mapping  (in  contrast  to  the  gen¬ 
eration  of  statistically  equivalent  networks)  ensures  that  subtle 
spatial  correlations  present  in  the  original  material  arc  retained  in 
the  network  structure. 

We  refer  to  the  network-generation  technique  presented  here  as 
a  grain-based  algorithm ,  whieh  refers  to  the  fact  that  the  first  step 
in  the  algorithm  (prior  to  network  generation)  is  the  characteriza¬ 
tion  of  the  underlying  granular  structure  in  the  material.  The  ben¬ 
efits  to  this  approach  arc  threefold.  First,  mapping  of  the  granular 
structure  is  less  sensitive  to  image  resolution  than,  for  instance,  the 
direct  mapping  of  a  pore  skeleton.  Second,  because  the  pores  in  a 
granular  material  arc  necessarily  defined  by  their  surrounding 
grains,  a  map  of  the  grain  structure  is  an  ideal  template  for  gen¬ 
erating  the  pore  network;  in  principle,  the  map  of  the  grain  struc¬ 
ture  should  be  essentially  independent  of  image  resolution,  which 
means  that  the  network  structure  will  correlate  more  strongly  to 
granular  structure  than  to  voxel  size.  Finally,  because  the  number 
of  grains  is  typically  orders-of-magnitude  smaller  than  the  number 
of  voxels,  the  grain-based  algorithms  are  faster  than  equivalent 
voxel-based  algorithms. 

In  this  work,  we  use  a  scries  of  computer-generated  sphere 
packings  to  validate  and  quantify  errors  in  the  grain-based  net- 
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TABLE  1— PARAMETERS  USED  TO  DEFINE  THE  PORE  NETWORK 


Variable  Association 

Variable  Name 

Variable  Type 

Dimension 

Network 

Domain  dimensions 

Veetor 

Length 

Pore 

Location 

Veetor 

Length 

Void  volume 

Sealar 

Length3 

Maximum  inscribed  radius 

Scalar 

Length 

Throat 

Intereonnectivity-periodieity 

Sealanveetor 

Cross-seetional  area 

Scalar 

Length2 

Maximum  inscribed  radius 

Sealar 

Length 

Surface  area 

Sealar 

Length2 

Hydraulie  conductivity 

Scalar 

Length3 

work-generation  algorithm.  The  algorithm  is  then  used  to  perform 
network  modeling  on  two  different  facies  from  a  well- 
characterized  outcrop  (Gani  and  Bhattacharya  2003),  which  illus¬ 
trates  the  potential  of  network  modeling  to  be  used  in  a  multiscale 
framework.  Few  details  of  the  grain-reconstruction  aigorithm  (the 
precursor  to  the  network-generation  algorithm)  are  given  in  this 
paper  because  they  can  be  found  elsewhere  (Thompson  et  al. 
2006).  It  was  originally  developed  for  the  analysis  of  marine  sands 
(Reed  et  al.  2005),  but  is  being  applied  to  other  materials.  This 
work  extends  the  grain-based  approach  to  relatively  low-porosity, 
cemented  sandstones,  which  pose  different  imaging  and  analy¬ 
sis  challenges  than  unconsolidated  sands;  however,  the  main  focus 
of  the  current  paper  is  the  generation  of  pore  network  models 
rather  than  grain  characterization.  Although  grain  structure  pro¬ 
vides  an  initial  template,  ultimately  the  parameters  in  the  pore 
network  are  computed  directly  from  the  voxel  image.  Hence  (as 
explained  next),  even  if  the  ievel  of  consolidation  in  a  sandstone 
presents  problems  for  the  grain-reconstruction  algorithm,  the  re¬ 
sulting  pore  network  will  stiii  be  a  quantitative,  one-to-one  map  of 
the  pore  structure. 

Background 

XMT.  XMT  provides  nondestructive  and  noninvasive  3D  images 
of  the  interior  of  objects  by  mapping  the  X-ray  absorption  through 
the  sample.  The  amount  of  absorption  depends  on  the  chemical 
composition  and  structure  of  the  material  and  the  X-ray  energy. 
XMT  is  based  on  the  reconstruction  of  the  cross-section  of  an 
object  from  its  projection  data,  which  is  obtained  by  passing  a 
series  of  rays  through  an  object  and  measuring  the  attenuation  of 
these  rays  using  detectors  placed  on  the  downbeam  side  of  the 
object.  Projections  are  obtained  by  measuring  the  X-ray  attenua¬ 
tion  coefficients  of  the  sample  at  different  angles  as  the  sample  is 
rotated  about  the  vertical  axis.  These  attenuation  values  are  rep 
resented  in  images  as  discrete  elements  (pixels  in  2D  images  and 
voxels  in  3D  images).  Synchrotron  radiation  has  several  advan¬ 
tages  over  traditional  X-ray  sources,  including  high  flux  Intensity 
(number  of  photons  per  second);  a  high  degree  of  collimation 
(source  divergence  leads  to  image  blur);  and  the  ability  to  tune  the 
photon  energy  to  a  single  energy  or  frequency  over  a  wide  range 
using  an  appropriate  monochromator,  which  can  be  used  to  make 
element-specific  measurements. 

Over  the  last  two  decades,  XMT  has  played  an  increasingly 
important  role  in  the  characterization  of  porous  media  flow  and 
transport.  Because  of  its  nondestructive  nature  and  increasingly 
high  resolution,  synchrotron-based  XMT  provides  the  high-quality 
datasets  necessary  to  capture  the  3D  microstructure  of  the  media 
(Liang  et  al.  2000a,  b;  Lindquist  et  al.  2000,  Al-Raoush  and  Will- 
son  2005b)  and  the  distribution  of  fluids  within  the  pore  space 
(Seright  et  al  2002;  Seright  et  al.  2003;  Al-Raoush  and  Wilison 
2005a).  The  ability  to  characterize  and  correlate  the  void  space 
microstructure  and  fluid  distributions  provides  data  to  improve  and 
validate  pore-scale  models. 

Network  Modeling.  Network  modeling  has  a  long  history  in  the 
oil  and  gas  industry,  beginning  with  the  landmark  paper  of  Fatt 


(1956).  For  an  extended  period,  network  modeling  techniques  em¬ 
ployed  lattice-based  networks,  usually  decorated  with  a  distribu¬ 
tion  of  pore-body  and/or  pore-throat  sizes.  These  lattice-based 
models  are  valuable  for  qualitative  studies  of  transport  in  inter¬ 
connected,  heterogeneous  structures.  However,  they  have  not 
proved  to  be  effective  for  quantitative  modeling  of  real  materials. 

Beginning  in  the  early  1990s,  new  techniques  were  developed 
for  quantitative  network  modeling  Bryant  et  al.  (1993b)  created 
physically  representative  network  models  from  the  highly  charac¬ 
terized  Finney  packing  (Finney  1970).  Oren  and  coworkers  created 
synthetic  (computer-generated)  sandstones  and  developed  a  tech¬ 
nique  for  extracting  networks  from  these  structures  (Bakke  and 
Oren  1997,  Oren  et  al.  1998).  Their  group  has  continued  to  de¬ 
velop  this  approach  and  the  resulting  networks  have  been  used  by 
a  number  of  other  investigators  (Patzek  2001;  Lopez  et  al.  2003; 
Valvatne  and  Blunt  2003).  Lindquist  et  al.  (1996)  worked  with  3D 
microtomography  images  and  computed  the  medial  axes  of  the 
pore  structure;  this  approach  was  then  extended  to  allow  for  direct 
generation  of  network  structures  (Sok  et  al.  2002).  Thompson  and 
Fogler  (1997)  applied  the  techniques  of  Bryant  et  al.  (1993b)  to 
computer- generated  packed  beds,  and  Al  Raoush  et  al.  (2003)  ex¬ 
tended  this  work  to  allow  for  network  structures  with  arbitrary 
connectivity.  Ioannidis  and  coworkers  have  used  simulated  anneal¬ 
ing  and  other  algorithms  to  produce  networks  that  conform  to  key 
statistics  in  real  materials  (Talukdar  et  al  2002,  Liang  et  al.  2000; 
Ioannidis  et  al.  1997). 

For  the  network  modeling  used  in  this  work,  we  borrow  the 
terminology  physically  representative  network  models  (Bryant  et  al. 
1993b)  to  describe  the  general  class  of  models,  and  we  note  two 
important  characteristics  of  these  structures.  First,  this  type  of 
network  is  a  one-to-one  mapping  of  a  specific  porous  material  of 
interest,  which  ensures  that  subtle  spatial  correlations  in  the  pore 
structure  are  retained.  Second,  the  networks  are  described  using 
rigorous  geometric  parameters,  which  ensures  that  the  pore  mor- 
phoiogy  is  not  compromised  despite  the  need  to  discretize  the  pore 
space.  This  latter  point  contrasts  with  techniques  in  which  the  pore 
structure  is  transformed  into  a  network  of  interconnected  capillar¬ 
ies  from  the  outset,  an  approach  that  has  been  shown  to  cause 
ambiguity  In  the  subsequent  modeling  of  flow  (Balhoff  and 
Thompson  2004). 

There  is  no  unique  or  correct  discretization  of  most  real  pore 
structures  (exceptions  being  simple  structures  such  as  cubic  pack¬ 
ings  of  spheres).  Likewise,  there  is  more  than  one  approach  that 
can  be  taken  to  describe  a  network  model.  In  this  work,  we  use  the 
set  of  parameters  shown  in  Table  1  to  describe  the  network  struc¬ 
ture.  In  addition,  the  network  remains  linked  to  the  original  data 
(whether  it  is  a  voxel  image  or  a  computer-generated  material). 
This  methodology  ensures  that  additional  characterization  could  be 
performed  if  warranted  by  a  particular  modeling  algorithm.  (It 
guarantees,  in  essence,  that  none  of  the  original  morphologic  data 
are  lost.) 

The  flow  modeling  itself  is  performed  by  imposing  conserva¬ 
tion  equation (s)  at  each  pore  in  the  network.  This  results  in  a  set  of 
linear  or  nonlinear  algebraic  equations,  depending  on  the  physics 
of  the  process  being  modeled.  A  description  of  flow-modeling 
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algorithms  is  beyond  the  scope  of  the  current  paper,  but  this  in¬ 
formation  can  be  found  in  many  other  papers  on  network  modeling 
(Bryant  et  al  1993b;  Bakke  and  0ren  1997;  Patzek  2001;  Lopez 
et  al.  2003;  Valvatne  and  Blunt  2003;  Thompson  and  Fogler  1997; 
Balhoff  and  Thompson  2004;  Bryant  et  al  1993a). 

Network  Generation  From  Voxel  Data.  Despite  the  advances 
described  in  the  previous  section,  what  remains  surprisingly  dif¬ 
ficult  is  the  generation  of  physically  representative  network  mod 
els  directly  from  microtomography  images  of  real  materials.  The 
difficulty  stems  from  the  distinct  differences  in  the  form  and  scale 
of  the  data  structures.  The  XMT  data  consist  of  hundreds  of  mil¬ 
lions  (sometimes  billions)  of  voxels  on  a  Cartesian  grid.  In  con¬ 
trast,  the  network  representation  of  the  same  porous  medium  is 
likely  to  be  described  using  a  much  smaller  data  set  (order  1 03— 1 0s 
pores),  which  are  interconnected  by  a  complex  structure  of  flow- 
paths  rather  than  aligned  with  any  specific  coordinate  system. 
Transforming  the  former  to  the  latter  is  a  nontrivial  procedure. 

Most  previous  work  in  this  area  employs  the  medial  axis  as  a 
basis  for  characterizing  the  pore  structure.  For  discretized  images 
such  as  those  from  XMT,  the  medial  axis  is  obtained  by  thinning 
the  original  pore  structure  until  a  one-voxel-thick  skeleton  re¬ 
mains,  or  by  computing  distance  transforms  in  the  void  space  (to 
the  void/solid  interface)  and  retaining  the  skeleton  of  local 
maxima.  Different  medial  axis  structures  are  obtained  depending 
on  the  method  used,  the  order  for  thinning,  the  type  of  distance 
transform  used,  and/or  the  rules  invoked  to  ensure  that  topology  is 
retained  (Lohman  1998).  Lindquist  et  al  (1996)  used  the  medial 
axes  of  3D  images  to  compute  statistical  parameters  associated 
with  the  void  space.  They  later  extended  their  algorithm,  using  the 
medial  axis  as  a  basis  for  obtaining  pore-  and  throat  size  distribu¬ 
tions,  which  in  turn  were  mapped  onto  a  network  structure 
(Lindquist  et  ai.  2000).  Specificaiiy,  the  mediai  axis  was  trimmed 
down  to  the  percolating  fraction  and  nodes  were  then  merged  (by 
comparing  the  distance  separating  neighboring  nodes  to  the  dis¬ 
tance  to  local  surfaces)  to  define  pore  locations.  Pore  throats  were 
found  by  dilating  the  medial  axis  until  the  dilated  cylinder  con¬ 
tacted  the  bounding  surfaces;  pore-throat  geometries  were  then 
obtained  by  triangulating  between  the  medial  axis  and  voxels 
along  the  perimeter  of  the  constriction.  Sok  et  el.  simulated  im¬ 
miscible  displacement  on  networks  generated  by  this  algorithm 
and  compared  them  to  immiscible  displacements  on  regular  lat¬ 
tices  having  identical  values  for  key  statistical  metrics  (Sok  et  al. 
2002).  Mirroring  what  Bryant  et  al  (1993a)  observed  for  single 
phase  flow,  they  note  that  multiphase  flow  behavior  is  not  repro¬ 
duced  correctly  on  statistically  equivalent  networks,  emphasizing 
the  need  for  algorithms  that  capture  true  pore  morphology.  Delerue 
et  al.  (1999)  applied  a  medial  axis  technique  to  a  3D  image  of  a 
resin  impregnated  soil,  and  then  defined  the  pore-size  distribution 
in  the  soil  by  measuring  the  maximum  inscribed  balls  for  all  voxels 
contained  on  the  skeleton.  Mercury  intrusion  was  simulated  di¬ 
rectly  on  the  voxel  map  (rather  than  a  network  representation). 
Delerue  and  Perrier  (2002)  describe  in  detail  the  various  compu¬ 
tational  elements  used  in  the  algorithm.  Silin  et  al  (2003)  em¬ 
ployed  a  similar  approach,  except  that  the  maximal  inscribed  bail 
was  found  for  each  void-phase  voxel  in  the  packing.  Though  com 
putationally  more  intensive,  this  approach  allows  the  pores  to  be 
found  independently  of  the  skeleton.  They  tested  their  algorithm 
using  computer-generated  sphere  packings  and  a  CT  image  of 
Fontainbleau  sandstone.  However,  pore  network  models  were  not 
created,  and  a  quantitative  assessment  of  the  pore  locations,  sizes, 
and  connectivity  was  not  made. 

Materials  and  Methods 

Sandstone  Samples.  The  samples  used  for  this  application  were 
taken  from  the  Wall  Creek  Member  of  the  Cretaceous  Frontier 
Formation,  Wyoming,  USA  as  part  of  an  integrated  geologic,  geo¬ 
physical,  and  engineering  study  (Gani  and  Bhattacharya  2003). 
Ten  wellbores  were  drilled  through  the  Wall  Creek  Member  and 
sample  plugs  were  selected. 

Sample  A  is  from  a  tidally-reworked  sandstone  facies  within  a 
delta  front  sandstone  body.  The  permeability  of  this  sample  is 


severely  reduced  by  calclte  concretions.  Concretions  locally  re¬ 
duce  the  permeability  by  several  orders  of  magnitude  in  approxi¬ 
mately  15  volume  %  of  the  tidal  sandstone  facies  (Lee  et  al.  2007) 
The  sediments  preserved  in  these  rocks  were  derived  from  uplifts 
to  the  north  and  west.  Lee  et  al.  (2007)  report  that  the  average 
composition  is  approximately  51%  quartz,  21%  rock  fragments, 
and  28%  feldspars,  with  calcite  as  the  dominant  diagenetic  cement. 
Feldspar  and  volcanic  components  have  been  highly  altered  by 
diagenesis.  In  uncemented  rocks  of  the  tidal  sandstone  facies,  po¬ 
rosity  averages  approximately  0.22  and  permeability  is  in  the  range 
of  tens  to  a  few  hundred  millidarcies,  with  a  mean  of  110-140  md 
(higher  for  fluvlal-dominated  deposits  of  the  delta  front). 

Sample  B  is  from  a  channel  sandstone  facies  within  the  same 
delta  front  sandstone  body  as  sample  A.  The  provenance  of  this 
sample  Is  thus  similar  to  Sample  A,  except  for  calcite  cement 
replacing  some  of  the  pore  space.  This  sample  is  interpreted  to  be 
from  a  relatively  mud-free  cross-stratified  channel  sandstone  fa¬ 
cies.  Depositional  considerations  would  imply  very  high  reservoir 
quality.  The  grain  size  is  on  average  larger,  and  sorting  is  better 
than  sample  A  However,  it  is  from  a  cemented  region  within  the 
sandstone,  and  its  reservoir  quality  is  therefore  relatively  low  The 
greater  abundance  of  calcite  concretions  in  some  channel  sand¬ 
stones  is  inferred  to  be  related  to  the  abundance  of  mud  chips  (or 
shale  clasts).  Although  mud  chips  would  decrease  initial  perme¬ 
ability,  their  inhibiting  effect  on  concretion  formation  actually  re¬ 
sults  In  higher,  post-diagenetic  reservoir  quality  for  channel  sand 
stones  with  mud  chips  (Lee  et  al.  2007). 

XMT  of  the  Sandstone.  Small  samples  of  the  sandstones  were 
impregnated  with  an  epoxy  resin  under  vacuum  and  then  cored  to 
a  length  of  25  mm  and  a  diameter  of  5  mm.  5-mm-long  sections  of 
the  cores  were  imaged  at  33.07  keV  energy  at  the  13-BMD  tomog¬ 
raphy  beamline,  operated  by  GeoSoiiEnviroCARS  (GSECARS)  at 
the  Advanced  Photon  Source.  Image  reconstruction  was  performed 
using  algorithms  developed  by  GSECARS  to  convert  CT  attenu¬ 
ation  data  to  3D  volumetric  data  The  resultant  3D  gray  scale  im 
ages  have  a  voxel  resolution  of  7.63  microns.  Segmentation,  the 
process  of  converting  a  gray-scale  Image  to  a  1  -bit  or  binary  image 
by  separating  the  images  into  two  populations  based  on  gray-scaie 
values,  was  performed  using  the  indicator  kriging  technique  (Oh 
and  Lindquist  1999)  that  is  part  of  the  3DMA  software  package.  In 
this  work,  the  pore  voxels  are  assigned  a  vaiue  of  0  and  the  grain 
voxels  are  assigned  a  value  of  1.  Fig.  1  contains  images  of  the 
binary  volume  files. 

The  combination  of  mineralogy  (i.e.,  higher  absorbing  ele¬ 
ments)  and  X-ray  energy  created  a  relatively  high  signai-to -noise 
ratio.  This  less-than  ideal  image  quaiity  made  it  difficult  to  com¬ 
pletely  resolve  the  solid  and  void  phases  and  to  remove  the  ring 
artifacts.  Since  the  time  that  these  experiments  were  performed, 
these  issues  have  been  addressed,  which  means  that  images  from 
current  experiments  are  generally  free  from  these  problems. 

Network  Construction 

Overview.  The  Input  for  the  algorithm  is  a  binary  volume  file 
describing  the  porous  material  of  interest:  voxel  labels  indicate 
whether  the  voxel  is  contained  in  the  void  phase  or  solid  phase  of 
the  material.  The  algorithm  operates  by  identifying  grains,  search¬ 
ing  for  pores  based  on  the  grain  locations,  and  then  creating  the 
interconnected  network  using  a  novel  restricted -bum  algorithm. 
The  significant  differences  between  this  algorithm  and  other  net¬ 
work  generation  techniques  include  the  following: 

1. The  first  step  in  the  algorithm  is  the  identification  of  the 
granular  structure,  a  process  that  is  less  sensitive  to  image  resolu¬ 
tion  than  techniques  such  as  skeletonization.  Consequently,  net¬ 
works  created  for  the  same  material  at  two  different  image  reso¬ 
lutions  have  very  similar  structure  and  properties. 

2.  The  algorithm  uses  the  granular  structure  as  a  template  to 
help  define  pores.  This  means  that  the  computationally  expensive 
search  for  pores  becomes  linked  to  the  number  of  grains  in  the 
image  rather  than  the  number  of  voxels,  and  the  computational 
penalty  for  increasing  image  resolution  is  less  severe  than  with 
other  algorithms. 
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Fig.  1 — Binary  volume  files  for  two  sandstone  samples:  (a)  tid¬ 
al  ly-reworked  sandstone  facies;  (b)  channel-sandstone  facies. 

3.  Because  the  grain  locations  provide  the  framework  for  de¬ 
termining  pore  structure,  the  algorithm  operates  without  ever  com¬ 
puting  a  skeleton  of  the  pore  space.  By  eliminating  this  step,  a 
number  of  problematic  issues  such  as  nonuniqueness,  dependence 
on  image  resolution,  and  the  formation  of  internal  loops  are  avoided. 

4  A  key  intermediate  step  between  the  original  binary  map  and 
the  finai  network  structure  assigns  an  integer  iabel  to  every  voxel 
(both  void  and  solid  phase)  to  mark  the  grain  or  the  pore  to  which 
that  voxel  belongs.  The  reason  for  emphasizing  this  intermediate 
step  is  that  it  makes  obtaining  statistical  information  and  construct¬ 
ing  the  network  straightforward  and  unambiguous. 

Grain  Reconstruction  Algorithm.  Step  1.  Dual-Phase  Burn .  A 
simultaneous  grain  phase  and  void-phase  bum  is  performed  (using 
the  terminology  of  Lindquist  et  al.  1 996).  Voxels  in  the  solid  phase 
are  labeled  with  positive  integers  denoting  the  bum  level  and  vox¬ 
els  in  the  void  phase  are  labeled  with  the  negative  of  the  burn  level. 
This  convention  is  not  necessaiy  if  the  bum  map  remains  coupled 
with  the  binaiy  material  array;  however,  the  use  of  opposing  signs 
allows  the  bum  numbers  to  be  written  over  the  initial  material 
array  without  losing  the  phase  information 


Step  2 .  Location  of  Extrema,  A  search  is  performed  to  find 
local  maxima  in  the  bum  assignments  These  local  extrema  are 
islands  of  one  or  more  voxels  that  are  surrounded  by  voxels  with 
lower  bum  numbers.  In  the  simplest  form  of  the  algorithm,  the 
local  maxima  are  taken  to  be  the  grain  centers.  However,  refining 
the  grain  locations  using  optimization  leads  to  better  results 
(Thompson  et  al.  2006).  Simultaneously,  the  local  minima  can  be 
found  (which  according  to  the  sign  convention  are  in  the  void 
phase)  In  practice,  this  step  is  skipped  to  reduce  computation  time 
because  the  void -phase  extrema  are  not  used  to  define  the  pore 
structure.  However,  the  minima  are  found  and  reported  for  com¬ 
pleteness  in  this  paper. 

Step  3.  Tessellation  of  the  Grain  Centers.  A  periodic  Delaunay 
tessellation  is  performed  using  the  locations  of  the  grain  centers 
(Thompson  2002).  The  purpose  of  the  tessellation  is  to  identify 
likely  pore  locations  based  on  the  granular  structure.  The  role  of 
the  tessellation  is  somewhat  subtle,  and  the  following  points  are 
worth  clarifying: 

1.  In  contrast  to  other  techniques  where  the  Delaunay  tessella¬ 
tion  is  used  to  define  pore  structure  (Bryant  et  al  1993b;  Thomp 
son  and  Fogler  1997;  Ai-Raoush  et  ai.  2003;  Bryant  et  ai.  i993a), 
it  does  not  influence  the  structure  of  the  pores  in  the  current  al 
gorithm.  In  fact,  the  only  restriction  that  it  Imposes  is  on  the  total 
number  of  pores:  the  algorithm  will  not  find  more  pores  than  the 
number  of  tetrahedrons  in  the  tessellation.  This  limitation  should 
not  be  of  practical  consequence  because  the  tendency  of  the  De¬ 
launay  tessellation  is  to  identify  too  many  pores  (by  splitting  single 
pores  into  multipie  tetrahedrons;  see  Al  Raoush  et  al.  (2003)]. 
Furthermore,  this  restriction  can  be  relaxed  if  necessary  at  the 
expense  of  increased  computation  time  (see  Item  3). 

2.  Use  of  the  Delaunay  tessellation  to  aid  in  the  identification  of 
pore  locations  provides  an  advantage  over  erosion/dilation  tech 
niques  because  the  number  of  pores  that  are  iocated  for  a  given 
porous  medium  is  relatively  insensitive  to  voxel  resolution  This 
issue  is  demonstrated  later  in  the  validation  section 

3.  Finally,  the  tessellation  is  a  valuable  but  not  essential  part  of 
the  algorithm  An  alternative  approach  would  be  to  perform  the 
optimizations  described  in  the  next  step  beginning  with  every 
void  phase  voxel.  However,  even  for  the  relatively  small  sand¬ 
stone  datasets  used  here,  this  approach  requires  solving  -5,000,000 
nonlinear  optimization  problems  (in  contrast  to  -25,000  optimiza¬ 
tions  when  the  tessellation  is  used  as  a  template) 

Step  4.  Locating  Pores.  In  the  Introduction,  we  note  that  the 
division  of  void  space  into  pores  is  somewhat  arbitrary.  In  this 
work,  we  use  the  distance  function  d(x,y,$,  whose  value  gives  the 
minimum  distance  from  the  point  [xty,z\  in  the  void  space  to  a 
grain  surface  (Luchnikov  et  al.  1999),  and  then  define  a  pore  as 
any  local  maxima  in  d.  In  practice,  this  definition  corresponds  to 
an  accepted  definition  for  pores:  locations  of  maximum-diameter 
spheres  that  can  be  inscribed  into  the  void  space,  and  that  are 
constrained  from  movement  by  the  surrounding  soiid  phase 
(Scheidegger  1974).  In  Step  4,  these  maximum-diameter  inscribed 
spheres  are  located  by  performing  repeated  optimization  proce¬ 
dures  (to  maximize  d)  using  the  tetrahedrons  as  seed  locations  to 
start  the  optimizations. 

The  optimizations  themselves  are  performed  using  a  modified 
Poweli’s  method  (Press  et  al  1992),  which  is  a  direction  set 
method  effective  for  situations  where  gradients  in  the  objec¬ 
tive  function  cannot  be  calculated  directiy.  In  essence,  the  pro¬ 
cedure  repeatedly  performs  ID  line  minimizations.  Various 
schemes  are  available  to  choose  the  minimization  directions,  most 
based  on  estimating  new  conjugate  directions  from  the  histoiy  of 
the  optimization. 

As  each  extremum  Is  found,  the  voxel  containing  its  x,y,z  lo¬ 
cation  is  marked  with  the  pore  number.  If  a  different  seed  has 
already  converged  to  this  same  voxel  (though  the  coordinate  lo¬ 
cation  would  rarely  be  exactly  the  same),  a  new  pore  is  not  added 
to  the  list.  Hence,  Step  4  ends  having  generated  a  list  of  N  pore 
locations  and  inscribed  radii  (i.e.,  maximum  d  values)  along 
with  N  void- phase  voxels  labeled  with  the  corresponding  pore 
number.  At  this  point,  the  pores  are  no  ionger  tied  to  their  seed 
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tetrahedrons  and  there  is  no  further  need  for  the  Deiaunay  tessel¬ 
lation  in  the  algorithm. 

Step  5.  Pore  Merging.  A  viable  option  is  to  use  all  pore  loca¬ 
tions  identified  in  Step  4,  and  proceed  with  constructing  the  net¬ 
work.  However,  in  reai  cases,  many  pores  overlap  with  one  another 
by  a  significant  amount,  the  overlaps  being  caused  by  one  of  two 
reasons,  (a)  the  locai  pore  geometry  causes  two  independent  ex¬ 
trema  to  be  in  proximity;  (b)  two  different  seeds  have  led  to  the 
same  extremum,  but  numerical  error  or  optimization  tolerances 
have  caused  the  computed  extrema  locations  to  differ  by  at  least 
one  voxel. 

In  either  of  these  cases,  there  is  good  reason  to  merge  two 
largely  overlapping  pores  into  one.  Various  merge  criteria  can  be 
devised.  In  the  current  algorithm,  pores  are  merged  only  if  one 
inscribed  sphere  encompasses  the  center  of  a  neighboring  in¬ 
scribed  sphere.  The  location  and  radius  of  the  larger  inscribed 
sphere  is  used  as  a  seed,  and  a  locai  optimization  is  performed  once 
again  to  verify  the  location  of  the  merged  pore  With  the  adjusted 
pore  locations  (owing  to  the  reoptimization  of  merged  pores),  we 
have  found  it  advisable  to  make  another  pass  to  check  whether 
additional  pores  should  be  merged,  and  indeed  to  continue  this 
merge— > re-optimization  procedure  iteratively  until  no  more  pores 
are  merged.  Step  5  ends  with  the  same  information  as  Step  4  (pore 
locations,  radii,  and  corresponding  voxel  assignments),  but  with 
fewer  and  more  spatially  independent  pores. 

Step  6.  Grain  and  Pore  Assembly.  Step  6  is  the  key  interme¬ 
diate  step  mentioned  in  the  previous  overview;  the  assignment  of 
all  voxels  in  the  Image  to  one  of  the  grains  or  pores  identified  in 
Steps  2  or  5.  respectively.  This  step  is  performed  using  a  novel 
restricted  bum  algorithm,  which  is  described  in  more  detail  in  the 
context  of  grain  reconstruction  elsewhere  (Thompson  et  al.  2006). 
In  short,  it  assembles  a  cluster  of  voxels  together  that  are  tied  to 
one  of  the  grain  or  pore  locations  respectively,  in  accordance  with 
the  local  geometry. 

A  summary  of  the  logic  for  this  restricted  bum  procedure  (Step 
6)  is  the  following: 

1.  Set  the  minimum  bum  level  to  the  largest  (absolute)  value 
found  in  the  pore  space. 

2.  Loop  through  all  voxels  in  the  domain. 

3.  If  the  current  voxel  borders  a  voxel  already  assigned  to  a 
particular  pore  and  its  absolute  burn  level  is  greater  than  or  equai 
to  the  current  minimum  burn  level,  assign  it  to  the  same  pore  as  the 
assigned  neighbor. 

4.  If  any  new  voxels  were  assigned  during  the  last  pass,  go  to 
Step  2. 

5.  If  no  new  voxels  were  assigned  during  the  last  pass,  reduce 
the  minimum  burn  level  by  1;  go  to  Step  2 

Step  7.  Pore  Morphology  and  Network  Construction.  Once 
Step  6  is  complete,  determining  morphologic  parameters  and  con¬ 
structing  the  network  is  a  straightforward  process.  The  totai  vol¬ 
ume  of  each  pore  is  obtained  by  summing  the  volume  of  all  voxels 
assigned  to  that  particular  pore.  The  inscribed  pore  radii  are  al¬ 
ready  known  from  the  Step  4  and  Step  5  computations. 

In  our  definition  of  network  structure,  pore  throats  have  no 
volume  but  rather  are  defined  by  the  faces  where  two  pore- 
eiements  come  into  contact.  Hence,  the  pore-throats,  like  the  pores, 
have  rigorous  geometric  parameters  associated  with  them.  The  unit 
normal  for  the  pore-throat  interface  is  found  by  averaging  the 
orientation  of  all  voxel  interfaces  associated  with  the  given  throat. 
The  totai  cross-sectional  area  of  the  throat  is  then  found  by  sum¬ 
ming  the  areas  of  voxel  faces  at  this  interface  projected  onto  this 
local  unit  normal.  Use  of  the  projected  area  prevents  overesti¬ 
mating  the  area  because  of  the  staircase-iike  surface  created  by 
the  voxels.  The  inscribed  area  of  the  throat  is  found  by  determin¬ 
ing  the  largest  inscribed  sphere  whose  center  is  located  on  the 
throat  interface. 

Grain  surface  area  Is  also  assigned  as  a  throat  parameter  be¬ 
cause  it  affects  permeability  as  well  as  phenomena  such  as  absorp¬ 
tion  and  chemical  reactions  on  grain  surfaces.  The  surface  area  is 
assigned  by  estimating  the  surface  area  for  each  surface  voxel,  and 
then  assigning  that  element  of  surface  to  its  closest  throat.  This 
approach  ensures  that  the  surface  area  is  conservative  (i.e.,  the  sum 


of  all  pore  throat  surface  areas  equals  the  total  surface  area  in  the 
packing),  thus  providing  a  good  theoretical  foundation  for  its  use 
in  modeling. 

The  last  parameter  that  should  be  mentioned  is  the  pore-throat 
hydraulic  conductance  of  each  throat,  which  is  necessary  for  com 
puting  permeability  and  performing  dynamic  flow  simulations  In 
general,  the  conductance  is  computed  using  some  combination  of 
the  previously  mentioned  parameters,  through  generation  of  an 
equivalent  capillary.  In  this  work,  we  use  essentially  the  same 
approach  used  In  Thompson  and  Fogler  (1997)  (without  the  FEM 
computations).  Further  work  includes  evaluating  the  best  way  to 
compute  throat  conductivities  for  various  types  of  materials. 

The  network  itself  is  defined  by  mapping  out  the  connectivity 
of  each  pore,  which  is  defined  by  the  list  of  neighboring  pores  that 
share  voxel-voxel  contacts.  There  are  no  limitations  on  the  struc¬ 
ture  of  the  networks  thus  obtained  (i.e.,  coordination  numbers, 
etc.),  and  the  resulting  network  files  have  exactly  the  same  format 
as  network  descriptions  from  computer-generated  media  (Balhoff 
and  Thompson  2004).  Once  the  network  is  constructed,  plots  of 
pore-size  distribution,  throat-size  distribution,  and  coordination 
number  can  be  made  directly  from  the  data  in  the  network  file, 
without  having  to  return  to  the  large  voxeiized  data  sets. 

Validation.  The  grain-reconstruction  algorithm  has  been  tested 
extensively  using  computer-generated  packings  of  spheres  and  cyl¬ 
inders  as  well  as  XMT  data  from  unconsolidated  particles  and 
sands.  Separately,  the  network-generation  algorithm  has  been 
tested  using  a  series  of  computer-generated  packings  of  spheres. 
This  approach  is  valuable  for  three  reasons:  First,  the  pore  struc¬ 
ture  in  the  material  is  known  exactly;  hence,  the  validity  of  the 
network  generation  process  can  be  assessed  in  a  quantitative  man¬ 
ner.  Second,  the  computer-generated  structures  can  be  discretized 
at  arbitrary  voxel  resolution,  which  provides  good  benchmarks  for 
the  accuracy  that  can  be  expected  with  data  from  real  materials. 
Third,  the  computer-generated  data  are  free  from  noise  and  arti¬ 
facts.  which  allows  validation  to  be  focused  solely  on  the  network 
construction  (rather  than  imaging  and  segmentation  issues). 

Three  sphere  packings  were  used  for  validation:  a  cubic  pack 
ing,  a  rhombohedral  packing,  and  a  random  packing.  The  two 
regular  packings  span  the  full  range  of  attainable  porosities  for 
monodisperse  spheres  and  the  pore  structure  is  known  exactly.  The 
random  packing  is  more  representative  of  real  materials.  No  single 
network  is  correct  for  the  random  case.  However,  the  modified 
Delaunay  Tessellation  (MDT)  algorithm  (Al-Raoush  etai.  2003)  is 
a  fairly  rigorous  method  for  extracting  the  pore  structure  and  can 
be  used  for  comparison  to  the  networks  generated  from  the  vox¬ 
eiized  images. 

Table  2  explains  the  notation  used  for  reporting  the  statistical 
results.  Selected  quantitative  data  are  presented  in  Table  3.  (Data 
for  the  computer  generated  packings  are  in  units  of  sphere  diam¬ 
eters;  data  for  the  sandstone  shown  later  are  in  microns.)  Note  that 
for  the  three  sphere  packings,  the  top  row  of  each  section  contains 
parameters  obtained  from  the  MDT  networks.  These  parameters 
quantify  the  error  in  the  voxel-based  networks.  For  the  cubic  and 
rhombohedral  packings,  the  MDT  values  agree  with  theoretical 
values  from  a  unit-cell  analysis  and  are  exact  within  the  numerical 
tolerances  set  in  the  MDT  algorithm.  As  mentioned  previously, 
there  are  no  “correct”  morphologic  values  for  the  random  packing. 
However,  the  MDT  algorithm  is  an  excellent  benchmark  because 
it  is  based  on  Delaunay  analysis  (Finney  1970;  Mellor  1989)  and 
allows  for  variable  pore  connectivities  (which  are  not  possible  with 
the  unmodified  Delaunay  tessellation). 

Identification  of  Pores.  Column  #£Tis  the  number  of  extrema 
in  the  void-phase  burn.  This  parameter  is  the  most  logical  choice 
for  determining  pore  locations  if  the  burn  information  were  to  be 
used  directly  However,  notice  the  very  strong  dependence  of  these 
values  on  image  resolution  (at  least  for  the  noncublc  structures). 
For  the  rhombohedral  packing,  the  number  of  void  phase  extrema 
varies  from  four  to  328  as  image  resolution  is  increased  over  an 
order-of-magnitude  range  (the  true  number  of  pores  is  192).  For 
the  random  packing,  number  of  bum  extrema  varies  between  63 
and  803  as  resolution  varies  from  10  VPD  to  125  VPD  (the  best 
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TABLE  2— EXPLANATION  OF  PARAMETERS 

USED  IN  TABLES 

Parameter 

Description 

A„ 

Cross  sectional  area  of  a  pore  throat  (ave) 

A* 

Surface  area  assigned  to  a  pore  throat  (ave) 

Dg 

Grain  diemeter  (ave) 

DP 

Pore  diameter  (inscribed)  (ave) 

Dpt 

Pore-throat  diameter  (inscribed)  (ave) 

tiEG 

Number  of  extrema  in  the  grain-phase  burn 

UEV 

Number  of  extrema  in  the  void-phase  burn 

K 

Permeability  (subscripts  indicate  flow  direction) 

ipT 

Length  of  pore  throat  (ave) 

MBG 

Maximum  burn  number  in  the  grain  phase 

MBV 

Maximum  (absolute)  bum  number  in  the  void 
phase 

Ng 

Number  of  grains 

Np 

Number  of  pores 

PLE 

Average  error  (%)  in  the  computed  pore  locations 

np, 

Number  of  pores  initiaiiy  found  from  the 
tetrahedron  seeds 

UPM 

Number  of  one-to-one  pore  matches 
(reconstructed  packing  vs.  original) 

tiPu 

Final  number  of  pores  after  the  iterative 
merging/reoptimization  process 

#  Tet 

Number  of  tetrahedrons  used  as  seeds  for  pore 
locations 

Wvirtm 

Minimum  number  of  voxel  faces  required  for  a 
pore  throat  to  be  formed 

#  Vox 

Number  of  voxels  in  the  data  set 

vp 

Volume  of  a  pore  (ave) 

VPD 

Voxels  per  particle  diameter 

VR 

Voxel  resolution 

Z 

Pore  coordination  number 

£ 

Other 

Porosity 

ave 

Arithmetic  average 

Sd 

Standard  deviation 

estimates  for  number  of  pores  in  this  packing  are  -400).  Ciearly, 
the  void-phase  burn  provides  a  poor  indicator  of  pore  location, 
which  is  related  to  the  problems  with  skeletonization  that  were 
mentioned  previously. 

In  the  current  algorithm,  the  grain  structure  is  used  as  a  tem¬ 
plate  for  determining  pore  locations.  The  initial  seeds  for  pore 
locations  are  based  on  a  Delaunay  tessellation  of  the  grain  loca¬ 
tions.  Table  3  lists  the  number  of  Deiaunay  tetrahedrons  (see  col 
umn  #7ef).  The  number  of  largest-inscribed-spheres  obtained  from 
these  seeds  is  listed  in  column  Note  that  the  numbers  in  both 
of  these  columns  are  relatively  insensitive  to  image  resolution 
since  they  are  tied  to  the  grain  structure  rather  than  the  voxel  structure. 

The  cubic  packing  has  significantly  fewer  inscribed  spheres 
(i.e.,  pores)  than  tetrahedron  seeds.  The  reason  for  this  difference 
is  that  assembling  a  cubic  pore  requires  at  least  five  tetrahedrons. 
If  the  optimization  routine  were  exact,  all  five  seeds  would  con¬ 
verge  on  the  same  central  pore  location.  In  reaiity,  more  than  one 
seed  will  usually  converge  on  the  same  central  voxel,  and  is  listed 
only  once.  The  fact  that  %Pf  decreases  with  increasing  resolution  is 
somewhat  counterintuitive  (because  higher  resolutions  provide 
many  more  voxels  in  which  the  optimizations  can  land).  However, 
the  reason  for  this  effect  is  the  increased  accuracy  of  the  optimi¬ 
zation  procedure  as  resolution  increases  (which  is  a  consequence 
of  more  precise  distance-to-surface  calculations). 

The  more  important  value  Is  the  final  number  of  pores  after 
merging  (#P,J.  For  the  cubic  packing,  this  value  is  exact  for  all 


resolutions.  For  the  rhombohedral  packing,  it  is  exact  for  resolu¬ 
tions  of  16.7  VPD  and  above.  For  the  random  packing,  an  exact 
number  of  pores  can  never  be  defined  unequivocally.  However,  the 
values  obtained  from  the  current  algorithm  and  the  MDT  algorithm 
agree  reasonably  well.  Results  from  the  random  packing  also  Il¬ 
lustrate  the  efficiency  of  the  grain-based  approach.  Examining 
the  number  of  tetrahedrons  (each  of  which  generates  a  seed  point 
to  search  for  a  pore)  shows  that  -620  nonlinear  optimization 
procedures  were  performed  to  find  the  initial  pore  locations, 
independent  of  image  resolution.  In  contrast,  -32  million  of  these 
same  computations  wouid  have  been  required  for  the  100  VPD 
image  if  each  void  phase  voxel  were  tested  to  find  its  maximum 
inscribed  sphere. 

Accuracy  of  Pore  Locations  and  Sizes.  The  pore  location  is 
defined  as  the  center  of  the  largest  sphere  that  can  be  inscribed  into 
a  given  void  space.  For  the  cubic  packing,  the  results  are  essen¬ 
tially  exact  In  cases  where  an  integer  VPD  value  is  used.  For 
noninteger  VPD  values,  the  voxels  cannot  be  divided  evenly  be¬ 
tween  pores  in  the  packing.  Consequently,  there  remain  slight 
distributions  in  pore  parameters,  even  at  high  VPD  values.  How¬ 
ever,  error  remains  fairly  iow  in  key  parameters  such  as  pore 
volume  (see  iy  and  pore  location  [Table  3,  column  PLE).  The 
rhombohedral  packing  is  a  much  more  rigorous  test  because  the 
pores  are  not  symmetric  with  respect  to  the  Cartesian  voxel  grid 
and  also  because  they  are  small.  The  error  in  pore  locations  is 
significant  at  low  resolutions  (where  pore  dimensions  can  be  as 
small  as  a  single  voxel),  but  it  appears  to  decrease  monotonically 
with  increasing  voxel  resolution. 

The  accuracy  for  the  random  packing  network  is  somewhat 
harder  to  assess,  mostly  because  the  number  of  pores  found  does 
not  agree  exactly  for  the  MDT  vs.  voxel-based  networks.  Column 
#PA/  in  Table  3  provides  statistics  for  the  fraction  of  pores  with  a 
one-to-one  match  with  the  MDT  network  (defined  when  the  center 
of  one  and  only  one  pore  from  the  voxel  network  lies  inside  an 
inscribed  pore  of  the  MDT  network).  Also  shown  is  the  average 
error  in  pore  location  for  these  one-to-one  matches. 

Pore  parameters  are  computed  with  increasing  accuracy  as  the 
voxel  resolution  increases.  This  trend  is  easy  to  confirm  by  exam 
ining  the  results  for  the  cubic  packing  The  same  trend  occurs  for 
the  other  two  packings,  but  Is  harder  to  assess  from  simple  statis¬ 
tics:  a  rhombohedral  packing  contains  two  distinct  types  of  pores 
(which  is  why  the  standard  deviation  remains  non-zero  for  both  the 
MDT  results  and  the  high-VPD  results)  and  die  random  packing 
has  a  distribution  of  pores. 

Coordination  Number.  Determining  the  pore  coordination 
number  (the  number  of  throats  emanating  from  a  pore  and  con¬ 
necting  to  other  pores)  remains  the  biggest  problem  for  this  algo¬ 
rithm  (as  well  as  other  network  generation  algorithms).  For  the 
cubic  packing,  the  coordination  number  is  exact  for  integer  VPD 
values,  again  because  of  the  fact  that  pore  boundaries  are  coinci¬ 
dent  with  voxel  boundaries.  For  resolutions  that  produce  noninte¬ 
ger  VPD  values,  the  average  coordination  numbers  are  slightly 
higher  and  the  maximum  coordination  numbers  are  significantly 
higher  than  they  should  be.  The  problem  is  illustrated  in  Fig.  2, 
which  shows  two  pores  assembled  from  a  cubic  packing  with  a 
noninteger  resolution  of  21.3  VPD  Although  they  are  constructed 
correctly  (within  the  limitations  of  the  voxel  resolution),  the  fact 
that  the  pore  boundaries  are  not  coincident  with  voxel  boundaries 
creates  a  contact  between  two  voxel  faces  at  a  position  that  should 
be  a  comer-comer  contact  in  reality  (and  therefore  should  not 
register  as  a  throat  connection).  The  numbers  from  the  noninteger 
VPD  cases  indicate  that  this  problem  occurs  less  often  than  once 
per  pore  on  average  (though  certain  pores  are  assigned  many  extra 
neighbors,  as  evidenced  by  the  high  values  for  the  maximum  co¬ 
ordination  numbers).  For  the  rhombohedral  packing,  the  problem 
is  similar,  though  no  resolution  generates  the  exact  coordination 
pattern  because  the  rhombohedral  pore  boundaries  never  lie  along 
voxel  boundaries. 

For  random  pore  structures,  the  problem  is  exacerbated.  The 
pore  geometries  are  more  irregular  (anomalies  are  not  limited  to 
comer-comer  contacts  as  with  the  cubic  packing)  and,  like  the 
rhombohedral  packing,  the  Cartesian  voxel  system  does  not  allow 
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Fig.  2— Non-physical  pore-throat  connection  in  a  cubic  packing 
of  spheres. 


for  perfect  breaks  in  the  pore  geometry.  Estimates  of  average  pore 
coordination  number  fall  into  the  mid-sixes  range.  These  values 
are  probably  too  high,  again  because  of  identification  of  voxel 
voxel  connections  that  should  not  register  as  pore-pore  connec¬ 
tions  (a  fact  also  evidenced  by  the  large  values  of  maximum  co¬ 
ordination  number). 

The  obvious  solution  to  this  problem  is  to  set  a  minimum  limit 
for  the  number  of  voxel -voxel  contacts  that  constitutes  a  pore 
throat.  For  the  cubic  packing,  requiring  pore  throats  to  be  com¬ 
posed  of  three  or  more  voxel  faces  led  to  perfect  coordination 
numbers  at  all  non-integer  resolutions  tested.  For  the  40  VPD 
random  packing,  Fig.  3  shows  a  histogram  of  the  number  of  voxel 
faces  found  at  a  throat  connection.  The  histogram  has  a  broad 
minimum  before  the  population  begins  to  increase  at  around  70 
voxel-faces  per  throat.  Interestingly,  the  theoretical  minimum  for 
throat  size  at  this  resolution  is  64  voxel-faces  per  throat,  which 
lends  credence  to  the  histogram.  These  two  factors  suggest  that 
connections  made  at  less-than  ~64  voxel-faces/throat  are  anoma¬ 
lies  and  should  be  discarded.  Table  4  shows  the  results  for  various 
minimum  voxel  limits  applied  to  the  40  VPD  random  sphere  pack¬ 
ing.  For  a  conservative  minimum  of  50  voxel-faces/throat,  the 
revised  throat  parameters  are  in  much  better  agreement  with  the 
MDT  values,  with  the  average  pore  coordination  number  reduced 
to  4.67.  A  more  aggressive  limit  of  63  voxel-faces  per  throat 
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Fig.  3 — Histogram  of  the  number  of  voxel-voxel  contacts  in  pore 
throats  in  a  random  packing  of  spheres. 


caused  the  other  parameters  to  deteriorate,  indicating  that  real 
throats  are  being  discarded  with  this  higher  limit. 

Unfortunately,  it  is  difficult  to  generalize  an  effective  rule  for 
how  to  limit  interconnectivity.  Creating  a  histogram  (e.g.,  Fig.  3) 
on  a  case-by-case  basis  is  a  sensible  and  fairly  easy  approach 
with  an  automated  algorithm.  However,  tomographic  data  from 
real  materials  that  we  have  tested  do  not  show  the  bimodal  distri¬ 
bution  found  in  Fig.  3  and  therefore  do  not  provide  a  strong  ra 
tionale  for  a  cutoff  value  We  suggest  the  following  approach.  For 
the  purposes  of  network  generation,  small  values  or  no  cutoff 
should  be  used  because  extra  throats  will  have  negligible  effect  on 
most  transport  processes  (because  of  their  very  small  size),  and 
because  one  risks  eliminating  what  the  tomography  has  identified 
as  “reaP  connections.  For  the  purposes  of  statistical  analysis,  the 
extra  throats  do  cause  a  problem  because  they  affect  the  calcula 
tion  of  coordination  number  and  average  throat  parameters.  In 
these  cases,  it  is  worth  the  extra  effort  to  assess  the  distribution  of 
voxel  contacts,  and  provide  a  cutoff  if  a  case  can  be  made  to  do  so. 
Finally,  the  use  of  geometric  averages  should  Improve  the  statis¬ 
tical  characterization  even  in  cases  where  small  throats  are  mis¬ 
takenly  identified. 


Results 

Networks  were  created  using  the  sample  A  and  sample  B  sand 
stone  images.  For  illustration  purposes,  ball-and-stick  depictions 
of  the  network  structures  are  shown  in  Fig.  4.  Comparing  these 
networks  to  the  Fig.  1  tomography  images  emphasizes  the  one-to- 


TABLE  4— CHANGES  IN  PORE-THROAT  PARAMETERS  BY  LIMITING  CONNECTIVITY 


Pore-Throat  Parameters 


Z 

Dpt 

Acs 

a5 

VFttin 

mm 

max 

ave 

ave 

Sd 

ave 

sd 

ave 

sd 

MDT 

2 

12 

4.60 

0.2627 

0  0846 

0.3548 

0.0963 

0.3369 

0.1233 

0 

t 

26 

6  60 

0,2083 

0.1015 

0.2484 

0.1479 

0.2392 

0.1448 

1 

t 

26 

6.60 

0.2083 

0.1015 

0.2484 

0.1479 

0.2392 

0.1448 

2 

1 

23 

6.08 

02215 

0.0950 

0.2673 

0  1386 

0.2607 

0,1475 

3 

t 

21 

5.87 

0.2269 

0  0921 

0.2757 

0.1338 

0.2697 

0.1477 

4 

1 

20 

5.71 

0.2310 

0.0898 

02820 

0.1299 

0.2770 

0.1472 

5 

t 

20 

5.71 

0.2310 

0.0898 

0.2820 

0.1299 

0.2770 

0.1472 

to 

1 

19 

5.24 

02422 

0  0842 

0.3019 

0.1164 

03018 

0.1414 

50 

t 

18 

4.67 

0.2523 

0.0814 

0.3251 

0.1008 

0.3386 

0.1255 

63 

0 

17 

4.57 

0.3113 

0  0960 

0  2981 

0.1536 

0.3108 

0.1792 

June  2008  SPE  Journai 


i  7  i 


Fig.  4 — Ball  arid-stick  schematic  illustration  of  network  models 
corresponding  to  the  two  different  sandstone  samples. 

one  mapping  that  is  achieved  using  a  physically  representative 
network  modeling  technique. 

For  the  sample  A  data  set,  two  different  network  reconstruc¬ 
tions  have  been  performed.  The  first  uses  the  grain-reconstruc¬ 
tion  algorithm.  For  comparison,  a  second  run  was  performed  in 
which  the  search  for  pores  was  performed  directly  on  the  void- 
phase  voxels,  without  use  of  the  grain  reconstruction  and  the  De¬ 
launay  tessellation. 

All  runs  were  performed  on  IBM  Power5  8-way  1 .90GHz  p575 
machines,  which  are  operated  by  LSU’s  high-performance  com¬ 
puting  center.  The  grain-based  algorithm  required  -1  hr  of  run 
time,  while  the  voxel  based  network  required  ~24  hrs  to  run. 

Morphologic  Parameters.  Fig.  5  contains  an  image  of  the  grain 
reconstruction  process  for  sample  A,  which  is  included  for  general 
interest.  Shading  is  assigned  randomly  to  grain  numbers  to  give  a 
visual  indication  of  the  discretization  into  individual  grains.  (The 
missing  piece  at  the  bottom  left  corner  is  an  artifact  caused  by  the 
graphics  software.) 

In  previous  tests  of  the  algorithm  on  tomography  images  of 
unconsolidated  sands,  two  problems  with  the  grain-reconstruction 


Fig.  5 — Grain  reconstruction  of  the  sandstone  sample  A. 


process  were  identified.  The  first  problem  is  single  grains  being 
“broken”  into  more  than  one  piece.  This  problem  is  associated  with 
noise  in  the  data  and/or  the  misidentification  of  grain  centers 
(which  in  turn  stems  from  either  unusual  grain  shapes  or  limita¬ 
tions  in  the  voxei  resoiution  and  resuiting  bum).  The  second  prob¬ 
lem  is  the  opposite  situation,  the  identification  of  a  single  large 
grain  that  should  probably  be  broken  into  more  than  one  piece 
(although  this  decision  is  subjective  in  consolidated  materials). 
Both  problems  appear  to  be  more  pronounced  for  the  current  sand¬ 
stone  data  compared  to  unconsoiidated  sands.  Solutions  to  these 
problems  are  being  implemented,  but  are  not  yet  validated,  and 
thus  were  not  used  in  the  current  work 

Fortunately,  the  main  problem  with  the  current  image  appears 
to  be  the  division  of  single  grains  into  more  than  one  piece,  which 
will  not  affect  the  network  generation  process  in  a  detrimental 
way;  an  increased  number  of  grains  will  simpiy  cause  the  algo¬ 
rithm  to  search  for  more  pores.  However,  the  final  pore  locations 
are  determined  by  the  optimization  procedure  and  reflect  local 
geometry  regardless  of  whether  the  local  grains  were  reconstructed 
in  a  reasonable  way.  Hence,  the  grain  based  algorithm  generates  a 
physically  representative  network  even  if  problems  occur  in  the 
grain  reconstruction  process. 

Table  5  contains  selected  parameters  for  the  two  different  net 
work  generation  techniques  used  on  Sample  A.  Fig.  6  shows  his¬ 
tograms  of  the  inscribed  pore-size,  inscribed  pore- throat- size,  and 
pore  coordination  number  for  the  grain- based  network.  Fig.  7  is  a 
graphical  representation  of  the  two  different  network  structures  in 
a  small  interior  section  of  the  sandstone. 

The  most  striking  difference  between  the  two  networks  is  the 
number  of  pores,  with  the  grain-based  network  having  many  fewer 
pores  and  thus  longer  pore  throats  spanning  the  gaps  between  the 
pores.  This  is  a  consequence  of  the  logic  in  the  grain-based  algo¬ 
rithm,  which  uses  the  Delaunay  tessellation  of  the  grains  as  seed 
locations  for  potential  pores.  In  contrast,  the  voxel- based  algorithm 
searches  for  pores  from  each  voxel  location,  which  results  in  the 
identification  of  many  more  local  extrema  (i.e.,  locations  of  maxi¬ 
mum  inscribed  spheres)  Although  the  difference  in  pore  numbers 
is  large,  it  is  still  not  possible  to  declare  one  network  more  correct 
than  the  other  because  discretization  of  the  pore  space  is  an  arbi¬ 
trary  process  (within  reason  of  course)  The  key  is  that  in  both 
cases  the  network  is  a  one-to-one  mapping  of  the  pore  structure, 
and  in  both  cases  the  pores  are  defined  using  the  same  rigorous 
definition:  a  local  maximum  in  the  distance  function  dfoy.z).  The 
most  effective  illustration  of  these  points  is  to  examine  the  two 
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TABLE  5 — QUANTITATIVE  RESULTS  FROM  NETWORK  MODELING 

OF  THE  SANDSTONE  IMAGES 

Sandstone  A 

Sandstone  B 

Grain-Based  Network 

Voxel-Based  Network 

Grain-Based  Network 

ng 

2,334 

2,185 

Dg  (,um) 

101 

100 

NP 

10,768 

65,574 

9,063 

DP  (/urn) 

34 

18 

39 

Dpri/um) 

31 

26 

38 

Lpt  {/urn) 

77 

35 

92  8 

Z 

3.21 

3.16 

3.15 

e 

0.186 

0.187 

0,160 

Kkx  (mD) 

410 

430 

1630 

Kyy  (mD) 

530 

570 

850 

(mD) 

320 

360 

631 

networks  side  by  side,  as  in  Fig.  7,  and  note  that,  despite  differ¬ 
ences  in  the  details,  the  trends  in  network  structure  are  the  same 
(e.g.  the  appearance  of  large  pores  at  the  same  locations;  gaps  in 
the  network  structure  at  the  same  locations,  and  so  on). 

The  other  notable  point  from  the  data  presented  here  is  the 
relatively  sinali  difference  between  the  pore-throat  diameters  and 
the  pore  diameters,  which  is  a  consequence  of  the  more  tube-like 
pore  structure  that  can  occur  in  consolidated  sands.  In  fact,  in  the 
voxel-based  network,  the  average  pore  size  is  smaller  than  the 
average  throat  size.  Although  this  is  counter  to  the  traditional  pore 
and  pore-throat  model,  it  is  a  consequence  of  the  high  density  of 
pores  in  the  second  network.  The  void  space  comprises  strings  of 
largely  overlapping  pores  rather  than  distinct  pores  connected  by 
long  throats;  the  pore-throat  diameter  is  simply  the  size  of  the 
channel  at  the  point  where  two  inscribed  spheres  overlap,  and  thus 
is  not  necessarily  smaller  than  the  adjacent  pores. 

Flow  Modeling.  The  networks  were  used  to  model  single-phase 
low-Reynolds-number  flow  of  a  Newtonian  fluid.  The  flow  mod¬ 
eling  was  performed  after  cropping  100  fx m  from  ail  sides  of  the 
networks  so  as  to  avoid  edge  effects  that  are  caused  by  the  bound¬ 
ary  of  the  volume  data. 

A  pressure  gradient  was  imposed  in  one  of  the  coordinate  di 
rections  and  the  resulting  volumetric  flowrate  was  computed.  Per¬ 
meability  was  then  calculated  using  Darcy’s  law  (see  Table  5). 
This  process  was  repeated  for  all  three  coordinate  directions  de¬ 
fined  by  the  imaging  voiume.  Note  that  KrI  does  not  imply  orien¬ 
tation  with  respect  to  a  bedding  plane;  because  the  sample  was 
taken  from  a  small  cutting,  no  attempt  was  made  to  align  the 
tomography  image  with  the  geologic  strata. 

The  ability  to  obtain  dimensional  permeability  values  is  a  con¬ 
sequence  of  using  physically  representative  networks.  That  is,  the 
network  is  a  map  of  a  specific-sized  region  of  the  sandstone,  which 
allows  for  the  computation  of  dimensional  volumetric  flowrates 
and  cross-sectional  areas  for  flow.  Additionally,  because  of  the 
one-to-one  mapping,  the  model  has  no  adjustable  or  scaling  pa¬ 
rameters.  [The  pore-throat  hydraulic  conductivities  are  determined 
using  the  formulas  originally  developed  for  sphere  packings 
(Thompson  and  Fogler  1997)  and  therefore  are  not  treated  as  ad¬ 
justable  parameters.] 

The  similarity  in  permeabilities  for  the  two  different  networks 
created  from  sandstone  A  may  be  surprising,  considering  the  dif¬ 
ferent  network  structures.  However,  this  fact  is  again  a  conse¬ 
quence  of  using  physically  representative  networks;  in  the  network 
with  fewer  pores  and  pore  throats,  the  throats  are  longer  and  there¬ 
fore  exhibit  larger  hydraulic  resistance.  Put  simply,  everything 
comes  out  in  the  wash. 

Fig.  8  is  a  grayscale  mapping  of  single-phase  flow  through  the 
same  region  of  sandstone  A  shown  in  Fig.  7b  (the  voxel  based 
network).  The  pressure  gradient  is  from  left  to  right  for  this  simu¬ 


lation,  and  brighter  shades  indicate  higher  flowrates.  This  graphic 
shows  clearly  the  heterogeneity  in  the  flow  pattern  that  results 
from  the  pore  morphology,  and  demonstrates  the  conventionai 
wisdom  that  relatively  few  pores  carry  a  large  fraction  of  the  fluid. 
This  graphic  is  also  an  effective  illustration  of  the  rationale  for 
using  network  modeling,  as  opposed  to  simpler  effective  medium 
or  bundle-of-tubes  models.  Although  these  simpler  models  may 
predict  correct  permeabilities  given  the  proper  pore  statistics,  they 
do  not  capture  the  flow  heterogeneity  that  contributes  strongly  to 
processes  such  as  solute  dispersion 

Discussion 

A  grain-reconstruction  algorithm  (originally  developed  for  uncon 
solidated  sands)  has  been  applied  to  consolidated  sandstones.  With 
additional  development  and  testing,  we  believe  that  this  algorithm 
can  be  used  as  a  tool  for  quantitative  analysis  of  the  granular 
structure  in  sandstones.  This  paper  uses  a  grain-reconstruction  pro¬ 
cess  as  a  template  for  generating  a  physically  representative  net 
work  model  of  the  pore  structure.  The  main  advantages  of  this  new 
approach  are  that  the  network  is  tied  to  a  correct  characteristic 
scale  for  the  materials  (i.e.,  the  grain  scale),  the  algorithm  is  more 
than  an  order  of  magnitude  faster  (compared  to  an  equivaient 
voxel-based  algorithm),  and  the  resulting  networks  are  less  depen¬ 
dent  on  voxel  resolution. 

Future  research  should  examine  how  different  network  struc¬ 
tures  (i.e.,  two  different  networks  that  are  intended  to  represent  the 
same  material)  affect  the  modeling  of  transport  phenomena.  To 
begin  studying  this  issue,  we  ran  the  current  algorithm  in  two 
ways;  as  a  voxel -based  algorithm  and  as  a  grain- based  algorithm 
These  two  approaches  resulted  in  dramatically  different  network 
structures;  the  most  obvious  difference  is  in  the  number  of  pores, 
but  this  in  turn  affects  most  other  parameters,  including  coordina¬ 
tion  number,  throat  length,  and  throat  conductivity.  Interestingly, 
despite  the  fact  that  the  voxel -based  network  contained  six  times 
the  number  of  pores,  the  average  permeability  difference  (for  the 
three  directions)  is  only  7.5%.  This  result  is  encouraging  because 
it  means  that  the  physically  representative  network  generation  pro¬ 
cess  (which  is  free  of  adjustable  parameters)  is  working  as  in¬ 
tended  with  respect  to  single-phase  permeability.  Nonetheless,  we 
expect  that  multiphase  simulations  will  be  more  sensitive  to  factors 
such  as  pore  coordination  number  because  current  multiphase 
models  rely  more  on  rule-based  algorithms  compared  to  the  single¬ 
phase  models.  This  is  a  topic  of  ongoing  research. 

We  conclude  by  commenting  on  the  issue  of  image  resolution 
From  our  experience  with  grain-based  modeling,  we  suggest  a 
minimum  of  10  VPD  in  order  to  extract  reasonably  good  statistics 
for  the  solid-phase  structure  in  unconsolidated  materials  (Thomp¬ 
son  et  al.  2006).  However,  this  value  (10  VPD)  may  need  to  be 
significantly  higher  for  pore  networks.  Consider  that  when  three 
spheres  are  placed  in  contact  to  form  a  pore  throat  (i.e.,  the  ge- 
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Fig.  6 — Morphology  of  the  pore  space  in  sandstone  A:  (a)  pore- 
size  distribution  (inscribed);  (b)  pore-coordination-number  dis¬ 
tribution;  (c)  pore-throat-size  distribution  (inscribed). 

ometry  in  a  rhombohedral  packing),  the  inscribed  diameter  of  the 
passage  is  equal  to  0.155Z)sphi;rc.  Hence,  a  10  VPD  image  will 
provide  only  1-2  voxel  resolution  at  a  pore  throat,  which  is  the 
crucial  area  for  flow  modeling.  This  issue  is  exacerbated  when  a 
particle  size  distribution  exists,  or  with  consolidation.  Table  3  from 
this  paper  provides  a  good  starting  point  for  assessing  the  errors  in 
pore  parameters  that  stem  from  voxel  resolution  issues. 
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